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ABSTRACT 

We present a scheme for the analytic computation of renormalization 
functions on the lattice, using a symbolic manipulation computer language. 
Our first nontrivial application is a new three-loop result for the topological 
susceptibility. 



1. Introduction 

The formulation of quantum field theories on the lattice has been mainly motivated by 
the need to study observables which are not amenable to a perturbative treatment. Yet, 
since the first days of the lattice, it became clear that perturbation theory could not be 
completely done away with; many quantities of physical interest measured on the lattice 
are connected to their continuum counterparts through renormalization functions which, in 
most cases, can only be calculated perturbatively. At a time when Monte Carlo numerical 
results are becoming increasingly accurate, higher order calculations of these functions, 
leading to non- negligible corrections, are necessary to achieve a matching precision. 

In the present paper we report on a scheme which we have developed for doing pertur- 
bative calculations on the lattice, using a symbolic computer language. Various schemes for 
doing similar calculations in the continuum exist since many years now, starting with Velt- 
man's Schoonschip; on the lattice, the lack of Lorentz invariance and the non-polynomial 
nature of the action introduce several additional complications, which we will point out 
below. We are currently working in formulating our computational scheme into a package 
for general use; in what follows we will limit ourselves to highlighting the essential points, 
deferring a detailed presentation of our algorithms to a future publication. 

As a first nontrivial application we also present the calculation to three loops of the 
additive renormalization (perturbative tail) of the topological susceptibility |l|-f§]- This 
operator has been studied for a number of years by several groups, using different meth- 
ods |§-|9[ , and is currently still under investigation, in particular in actions with dynamical 
fermions [|10|-[lT| and around finite temperature phase transitions ; both the presence of 



a phase transition, and the need to test further for agreement among the methods adopted, 
call for more precise determinations of this operator. 

A final introductory remark is in order here: In dealing with perturbation theory, one 
must bear in mind some well-founded caveats, stemming from the asymptotic, non-Borel 
summable nature of the perturbative series. As an example, the task of subtracting 
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additive renormalizations (mixing with lower dimensional operators) from a Monte Carlo 
signal seems rather problematic in principle. While general (non-) feasibility proofs are 
lacking, there do exist demonstrations, at least in 2-dimensional vector models, that some 



consistently defined and unambiguously separated from the physical signal. At any rate, 
it goes without saying that consistency checks are very important in these calculations, to 
ascertain that numerical results do show the expected theoretical behaviour; fortunately, 
this has been the case with most observables considered so far. 

2. Lattice Perturbation Theory by Computer 

The tasks one must carry out in doing lattice perturbation theory on a computer (or 
otherwise) are, in a nutshell: 
a) Computing the vertices 

P) Generating all relevant diagrams (with correct weights) 
7) Performing the contractions for each diagram 

5) Extracting powers of external momenta from the resulting n-point function 

e) Producing numerical code for loop integrations. 

These tasks are independent of one another; in particular, one may choose to perform 
only some of them symbolically and the rest by hand. Although many of the issues involved 
are a standard part of lattice perturbation theory, we will highlight them in a way that 
points out to their algorithmic resolution. We will draw examples from the topological 
susceptibility, defined as: 



of these problems can be circumvented |TB[ ; thus, for example, a perturbative tail can be 





Q(x) being the topological charge density: 





Using a lattice version Ql(x) of Q(x): 
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Q(x) = lim-i-^Q L , (2.3) 
X can be obtained from Monte Carlo data of Ql through: 

pj^ £<0 \Q L (n)Q L (0) | 0)m.c. = X a 4 (^(/3)) 2 + (£^F^ u )a 4 b((3) + d{(3) (2.4) 

n 

The first nontrivial loop calculations of Z((3), b(/3), d{j3) were done in Refs. 0-0], PH] . In 
this work we calculate the 3-loop coefficient c?4 of c?(/3) =c?3 / /? 3 +ti4/ /3 4 H — • . 

Let us go briefly through the chief points in the above tasks. 

a) Since both the action and most operators on the lattice are written in terms of link 
variables, which are exponentials of the gauge potential, they contain vertices with an arbi- 
trarily large number of gluon lines. The size of these vertices grows in principle as ( n+ ^ +1 ) , 
where n is the number of gluon lines and / is the number of links in the corresponding 
operator; to give a rough idea, already at sixth order the vertex for the topological charge 
occupies several dozens of output pages. Thus, any algorithm for generating vertices must 
take great care to keep them compact. 

Generically, an n-point vertex can be written in the form (for the sake of simplicity, 
we shall omit throughout this presentation our treatment of ghosts and fermions; these 
present some further complications, but no real stumbling blocks): 

J d% . . . d±k n ^ (f[ Aa jM e~ ik ^ 2 j V^ih, ...,k n ) (2.5) 

a\,—,ari 

A^ik) is the gluon field with momentum k and Lorentz (color) index \x (a). The phases 
exp(— ik-^,/2) stands also for the unit 4- vector in the /i direction) are absent when the 
fields A^ik) are taken to reside on the center of the link; otherwise, they are explicitly 
pulled out of the vertex to ensure that Hermitian operators lead to vertices with Vs 
which are real (aside from an overall prefactor). Since such phases are carried by all 
vertices, including the propagator, they cancel out upon contraction in any gauge; thus all 
dependence on the location of each field within a link disappears. 

We obtain V either by a straightforward expansion of the exponentials or by iterative 
use of the Baker-Campbeil-Haussdorff (BCH) formula 

e A e B = e A+B+l\A,B\+- (2 .6) 
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which results in terms with more definite symmetry properties. Due care must be taken to 
assign different dummy indices and momenta to multiple powers of the fields; this is only 
one of many aspects of the computation which are trivial by hand, but not in an automatic 
evaluation on a computer. We put V in the form 

V = d (Li(£?i + E 2 + ...)+ L 2 (E 3 + E 4 + ...) + .. .) 

(2.7) 

+ C 2 (L 3 (E 5 + E 6 + ...) + L A {E 7 + E 8 + ...) + ...) + .. . 
where Cj are 'colour structures', Lj are 'Lorentz structures' and E{ are monomials in 
trigonometric functions of momentum components. As an example, for the 3-point vertex 
of Ql, V takes the form 

^(Jt b fc 2 , fc 3 ) = -32 % g 3 e ix < kl+k2+k ^ /™a 3 . (2.8) 

(^2 £ WM2M3Pi I 2 cos(fci 7/1/2) cos(/c r/ u 2 /2) cos(kyfj, 3 /2) cos(/c 2 -/U3/2) cos(/c 3 -^ 2 /2) sin(/c r pi) 
pi 

+4/3 cos(fci-^i/2) cos(/c 2 - / u 2 /2) cos(/c3-//3/2) sin(fei-/0i/2) sin(/c 2 -pi/2) sin(/c 3 -pi/2)] 

_ X] ^iM2 £ wM3pip 2 [ cos ( /c 3-/U3/2) cos(/c 3 -p 2 /2) sin((/c 2 -/ci)-p 2 /2) sin(/c 3 -^i/2) sin(/c 3 -pi)]} 

P1P2 

Beyond 3-point vertices more colour structures can arise, for example: ^ c f aia2C f a 3 a ic 
anc l ^2 c d aia2C d a3aiC ; even though, in principle, a single structure, tr(T ai T a2 ■ ■ -T a ™), would 
suffice for any vertex (T a is a generator of the gauge group), it is preferable to use more 
symmetric structures for the sake of compactness. Lorentz structures proliferate on the 
lattice due to lack of rotational invariance; they also require use of 'internal' Lorentz indices 
(denoted pi) which are summed over. 

At this stage it is crucial to exploit the fact that all vertices are completely symmetric 
under interchange of external lines. We use this symmetry to compactify the corresponding 
expressions for V in three steps: First, reduce colour structures to a minimum, e.g. put 

fa-ta^c fa^aic - n ^ e form J2 C f ai<l2C / a 3 a 4 c ; second, use the residual symmetry of each 
color structure to reduce all accompanying Lorentz structures to a minimum; and third, 
for each color-times-Lorentz structure use its residual symmetry to reduce the number of 
accompanying monomials to a minimum. 

Some other aspects of the construction of vertices are, in brief: Using up the symmetry 
under exchange of internal indices (pi) for compactness (this becomes more subtle when 
internal momenta are also present, as is the case with the effective vertices of Fig. 1 (a,b), 
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which are very convenient constructs); defining new 'tensors', such as 5 y 



not summed over ^2), together with their lists of tensorial properties; establishing a stan- 
dard, 'canonical' form for the trigonometric monomials in order to reduce their number 
to a minimum. This last issue is rather nontrivial and still lacks a satisfactory resolution; 
the point is that the most immediate candidate prescriptions for a canonical form, such 
as using monomials with only one trigonometric function for every direction or using only 
ki/2 as arguments of these functions, have precisely the opposite effect of what is desired, 
leading to unmanageably large expressions. 

Since, for any given operator, vertices need be constructed only once in the beginning 
and then stored for subsequent use, considerations of speed are rather marginal here; they 
become far more pressing in what follows. Considerations of RAM usage are the main 
concern, since they determine the feasibility of this step of the computation on a given 
computer. 

P) The algorithmic generation of diagrams, together with their numerical weights, is 
the task which most resembles that of the continuum, the only difference stemming from 
the plethora of lattice vertices. For this reason, we shall not dwell on our approach, noting 
also that among the five tasks on the outset this is the only one still feasible by hand (given 
that calculations allowed by present computer capabilities can hardly reach 5 loops). 

As for the numerical weight of any given diagram, we can readily compute it from the 
formula: 



Here, w ex p is the product of (— l) k /k\ for each group of k > 1 identical vertices coming from 
the exponential of the action. The index i runs over all vertices in the diagram; the ith 
vertex has a total of n{ legs, of which ei remain external, b„ (even) get contracted among 
themselves, b{j get contracted against legs of the jth vertex. Finally, n g is the number 
of vertices of type g in the diagram and is the cardinality of that subgroup of the 
permutation group of all identical vertices which leaves b%j invariant (acting simultaneously 
on its rows and columns). In fact, #5 is the only quantity in W not given by a closed 
formula; however, it is a trivial matter to generate it numerically from the 'incidence 
matrix' bin. 



'exp 




U 9 (n g l) 

#5 



(2.9) 
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The diagrams contributing to x at 3 loops are shown in Fig. 2. Absent from this list 
are circa 140 diagrams involving the 2-point vertex of Ql as well as the effective vertex of 
Fig. lc, since both these vertices vanish in the case of zero external momentum (provided 
no infrared divergences are generated in this case, which we explicitly check). 

7) From an algorithmic point of view, the contraction for a given diagram entails the 
following: 

i) For each vertex involved, look up the corresponding expression for V and symmetrize 
it partially, according to the incidence matrix; thus, if two vertices are connected by n 
lines (bij=n), only one of the two need be symmetrized with respect to those lines. In 
a diagram like that of Fig. 2k, this means a potential saving of a factor of 4! in memory 
for intermediate expressions. Similar considerations apply for ej > 2 and bu > 4. 

ii) Form the product of all partially symmetrized vertices, renaming all indices (and mo- 
menta) as follows: Indices assigned to contracted legs become internal (fii^pii, aj-^Cji, 
kk^Pk'), and both internal and external indices are placed in ascending order, so as 
to make sure that their names remain distinct. 

iii) For each element of fry (i>j) consider in pairs the first bij available powers of the gauge 

the propagator: 

■ $(Pi+Pj) S CiCj {pf 8 PiPj - (1 - a)p iiPi pj iPj ) (2.10) 

( ft* = e'** - 1, P? = &'A> 

p=l 

Similarly for self-contractions {pa > 2). 

iv) Simplify color structures. Using the identity: 

1 1 

T i a j T kl = 2 5il5 *3 ~ 2N 5ij5kl C 2 ' 11 ' 
(valid for SU(N)), all internal color indices are completely eliminated. The algo- 
rithms which implement this simplification are identical to the ones used in the con- 



field from the i th and the ] th vertex (say, A p \(pi) and A^-ipj)) and substitute them with 
(2tt) 4 . . 1 

—^o{pi+pj)d CiCj ^2 "/',/', - > 1 - '"/''-/<;/'./*/<, 



tinuum 15 



Eliminate all Kronecker S's involving internal Lorentz indices. Doing so requires a 
judicious partial expansion of the expression (which is a product of large sums), to 
avoid drastic increases in memory. 
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vi) Compactify the expression using the symmetry under permutations of the names of in- 
ternal momenta and Lorentz indices. Allowed permutations for momenta are those 
consistent with the topology of the corresponding diagram, that is, permutations 
which leave invariant the momentum conservation delta functions 5(p\+ ■ ■ ■ +p n ) (or 
exp ix-{p\-\ — • +p n )) at each vertex; a table of these permutations is constructed right 
from the beginning. A conceptually easy algorithm would now generate all permuted 
versions of each subexpression and then select the first version in some order (e.g. lex- 
icographic). The problem is that both intermediate memory and execution time will 
grow factorially with the number of indices; since we often encounter up to ten indices, 
already at three loops, it is clear that such an algorithm will not do. At the price of a 
rather complicated source code, we have come up with an ordering algorithm which is 
(quasi-) polynomial in nature. 

Other considerations made here are ordering momenta simultaneously with indices, 
and casting our expressions in a form involving at most four Lorentz indices (given 
that the theory is defined in four dimensions). 

vii) Finally, trigonometric simplifications are systematically performed throughout, in order 
to put all terms in a canonical form, which then makes identifications or cancellations 
automatic. 



8) Extracting the analytic, exact momentum dependence of an n-point function, in the 
limit a—>0, is one of the most complicated tasks, both conceptually and algorithmically. 
This task does not enter the calculation which we present here; its elaboration is still in 
progress, and will be essential for the calculation of multiplicative renormalizations. 

Even in continuum regularizations this problem is only completely resolved at one 
loop [|T6 ] , ||17|| . To arbitrary loops, no systematic analysis of n-point functions exists. On 
the lattice, the first step is to decompose a given expression (to be integrated over internal 
momenta) in terms of a limited set of potentially divergent integrands, plus other terms 
which can be evaluated by setting the external momenta directly to zero. A possible basis 
for this set is: 



0<a lfl < < ax (2.12) 
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This decomposition poses no conceptual difficulty, but can cause disproportionate increases 
in memory, unless it is carefully implemented. One must now integrate the above set, ex- 
pressing the result in terms of standard functions (logarithms, Spence functions, etc.) and 
numerical constants characteristic of the lattice. At one loop this has been done system- 
atically ]r3|] , using a dimensional regularization technique JOJ . At higher loops, not only 
do these integrals become quite complicated, but their number also grows significantly. We 
are presently developing algorithms for carrying out the integration of the basis functions 



2.12)| symbolically. 



e) Having arrived this far, the only remaining task is the numerical evaluation of loop 
integrals with no dependence on external momenta. We do this both for finite and infinite 
lattices. 

On a finite lattice, loop integrals become nested multiple sums (with due attention paid 
to propagator zero modes). A mere conversion of the integrand to Fortran or C syntax is 
almost immediate. However, to produce optimized code one must take into account the 
following factors: 

i) Under certain changes of variables ({pi/j.^Pw, Vi}, {pi^— »— Pi^, Vz}, etc.) the integrand 
stays invariant (or can be rendered invariant at a small expense in size). It thus suffices 
to integrate over a small hypertriangular region of the original domain {— n < pi^ < 
7r, Vi, //}. An added complication for finite lattices is that the boundaries of this region 
are sets of nonzero measure. 

ii) Most diagrams contain two or more loops with no propagator line in common. (Among 
three loop diagrams the 'Mercedes', Fig. 2f, is the only exception.) Integration over the 
corresponding loop momenta need not be nested, but can be done independently, since 
all denominators (propagators) can be factorized into terms containing at most one such 
momentum. In order to factorize numerators as well, one must expand trigonometric 
functions containing more than one of the above loop momenta (again, our algorithms 
are written with an eye on keeping expansions to a minimum). The computational 
load for the resulting code is comparable to that of lower loop integrals. 

iii) The trigonometric functions comprising each monomial in the integrand typically de- 
pend only on a very small subset of the integration variables and can thus be pulled 
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out of most nested integrals. Further, since such factors are shared among many mono- 
mials, one can organize them in an (inverse!) tree, to avoid redundant integrations. 
We have incorporated all the above considerations in an algorithm which takes an 
integrand as input and produces optimized Fortran code for its integration. For lattice 
sizes of interest (~16 4 ), this optimization results in a gain in execution time of a factor of 
10 7 ! 

On infinite lattices, a drastic optimization is achieved by putting all propagators in the 
Schwinger representation: 

= / e-^dai (2.13) 
Pi Jo 

In this representation, integrations over different spatial directions factorize, so that their 
effective number is reduced by 3L—P (L: # of loops, P: # of propagators). At least one 
of the remaining integrations can be done analytically in terms of Bessel functions, leaving 
fewer integrals to be done numerically and a less singular integrand. We illustrate this 
with a very simple example: 

— 7^1 = / da * das e~ 8 ^ +a ^ $ 4 («i, a 2 , a 3 ) (2.14) 

-if p 2 q 2 p+q ( 27r ) ( 27r ) Jo 

f +7T dpdq 

$(ai,«2,«3)= / 77T -To exp(2ai cosp + 2ct 2 cos q + 2a^ cos(p+q)) 
/_„ ( 27T) 

r+ ,\ ' . (2-15) 

= / /e 2Q2COSp / (2 X /a2+«3+2«i«3Cosp) 

(Jo is the modified Bessel function.) The resulting expressions are amenable to high- 
precision Gauss-Legendre type integration. 

We also compare results for the infinite lattice to an extrapolation of finite lattice (L 4 ) 
results, of the form: 

result (L) = A + — (2.16) 

In addition to A and B, the exponent n may also vary for different diagrams (and one may 
also expect logarithms of L, cf. [Ref. 20| ). The discrepancy between this extrapolation and 
the infinite lattice result is typically a fraction of one per mille. 

In Table I we present the results for each diagram contributing to d^, for lattices of 
different size. Adding up all the contributions, we obtain for d^ (on an infinite lattice): 

d 4 = iV 4 (N 2 - 1) (1.735iV 2 - 10.82 + 73.83/iV 2 ) 10~ 7 (2.17) 



10 



Our programs are implemented in the computer language Mathematica. All numerical 
computations are performed by separate Fortran programs generated by our Mathematica 
routines. The generation of all vertices of Ql with up to six legs (needed for the two-loop 
calculation of Z(/3)) requires approximately 30 hours on a SUN Sparc Station 2 with 32 
Mbytes of RAM. The computation has to be split into ~ 200 independently computed 
contributions (summed only at the very end) in order to fit intermediate results into 
available RAM. 



3. Results and Conclusions 

The value obtained for allows one in principle to extract x with greater precision. With 
this value and those previously obtained for d% we have performed a series of best fits of 
Eq. (2.4)| to Monte Carlo data for SU(2) and SU(3). As in Refs. [0, | we have neglected 



the mixing with (FF). The values obtained for Xu (the non-renormalized topological 
susceptibility, Xu=X Z 2 (P)) are shown in Table II. They are in good agreement with those 
of Refs. 0, |§ (q.v. for details). At this stage, increased statistics of Monte Carlo 
data would be quite welcome on three counts: Improving the estimate of x, assessing 
the importance of the (FF) mixing, and checking for nonperturbative contributions in d n 
(such contributions are possible only for sufficiently high n, cf. |Ref. 2 If) . 

In conclusion, the calculational scheme which we have developed allows us to perform 
lattice perturbative calculations automatically, with very little 'human intervention'. Our 
aim is to be able to repeat the computation for different lattice operators without further 
programming, and this has not yet been completely achieved. The greatest difficulty that 
had to be overcome was the existing constraints on computer time and memory, which ne- 
cessitate devising polynomial-type algorithms and optimizing every aspect of this scheme. 
One major task still left open is the algorithmic extraction of external momentum depen- 
dence. Our first original application was the evaluation to 3 loops of the perturbative tail of 
X- We have also obtained three loop results for the gluonic condensate, and report them in 
a forthcoming publication p2[ . Repeating these calculations in the presence of dynamical 



fermions is relatively straightforward: Only a few additional diagrams appear, requiring no 
further computational resources. The calculation of multiplicative renormalizations within 
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this scheme, as well as a technical description of our algorithms, are postponed to a future 
publication. 
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Table I 



We list here the contribution to of individual diagrams, shown in Fig. 2, in the 
Feynman gauge. We use an L 4 lattice and gauge group SU(N). Each entry must be 
multiplied by N 6 (N 2 - 1) xl(T 7 . 



Fig. 


L = 3 


L = 8 


L = 16 


L = oo 


a 


31.23-33.75/iV 2 


37.79-41. 25/iV 2 


37.79-41. 38/iV 2 


37.77-41.40/iV 2 


b 


-2.137 


-2.981 


-3.081 


-3.114 


c+d 


5.210 


10.07 


10.78 


11.02 


e 


-4.328 


-5.964 


-6.162 


-6.228 


f 


.3955 


.7233 


.7429 


.7457 


g 


.0882 


.3107 


.3100 


.3099 


h 


1.988 


1.063 


.8017 


.7082 


i+j 





9.242 xl0~ 6 


9.112 xl0~ 6 


9.109 xl0~ 6 


k 


12.36-18.85/iV 2 


17.06-24.49/iV 2 


17.56-24.60/iV 2 


17.73-24.61/A^ 2 




+56.56/iV 4 


+73.48/iV 4 


+73.81/iV 4 


+73.83/iV 4 


1 


-44.09+44.99/iV 2 


-56.37+55. 00/iV 2 


-57.02+55. 17/A^ 2 


-57.20+55. 18/A^ 2 



Table II 

We list the values of Xu/Aq CD (for gauge groups SU(2) and SU(3), as obtained from 
Eq. (2.4)1 an d Monte Carlo data 0, || through a series of fits, in which: a) d± was 
an additional parameter to be fitted, or: b) the exact value of was taken from our 
calculation, fitting instead d§. 





SU(2) 


SU(3) 


a 


2.35(22) 10 4 


2.58(64) 10 5 


b 


2.07(23) 10 4 


2.70(69) 10 5 
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a) 



b) 



c) 



Fig. 1 Effective vertices. Solid (dashed) lines represent 
gluons (ghosts). Bullets are insertions of Q. 





i) 






k) 



1) 



Fig. 2 Diagrams contributing to . A cross 
stands for the integration measure. 
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